{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "28dd791a-e84a-441d-9952-8a3dc37f5d17",
   "metadata": {},
   "source": [
    "<img src=\"../common/rfsoc_book_banner.jpg\" alt=\"University of Strathclyde\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbdfad41-c485-492a-a99f-25c2dca975aa",
   "metadata": {},
   "source": [
    "# Notebook Set F\n",
    "\n",
    "---\n",
    "\n",
    "## 01 - Introduction to Frequency Planning\n",
    "In this notebook, we investigate spurs resulting from the data conversion process and how they can be calculated. We first discuss how non-linear responses can cause harmonics, using an example of a non-linear amplifier, then go on to look at interleaving ADCs and how mismatches can result in a variety of spurs across the spectrum. We discuss the effects of aliasing and how it affects frequency planning for ADCs due to spurs folding into the Nyquist band. We then go on to simulate a non-linear interleaved ADC and compare the practical results with the theory. Finally, we touch on designing filters to sufficiently suppress spurs close to the signal of interest.\n",
    "\n",
    "## Table of Contents\n",
    "\n",
    "* [1. Introduction](#introduction)\n",
    "* [2. Effects of Non-Linearity](#non_lin)\n",
    "* [3. ADC Interleaving Spurs](#il_spurs)\n",
    "    * [3.1. DC Offset Interleaving Spurs](#dc_offset)\n",
    "    * [3.2. Gain/Time Interleaving Spurs](#gain_time)\n",
    "    * [3.3. Harmonic Interleaving Spurs](#harmonic)\n",
    "* [4. Aliasing](#alias)\n",
    "* [5. Simulating a Non-Linear Interleaved ADC](#sim)\n",
    "* [6. Filter Design](#filt_design)\n",
    "* [7. Conclusion](#conclusion)\n",
    "    \n",
    "## Revision\n",
    "* **v1.0** | 05/12/22 | *First Revision*\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cee53fce-a6ad-4db4-a944-8b504ccec30d",
   "metadata": {},
   "source": [
    "## 1. Introduction <a class=\"anchor\" id=\"introduction\"></a>\n",
    "Frequency planning is a radio design technique used to avoid spurious emissions that may interfere with a signal of interest. It is based on the idea that the frequency content of many spurs are deterministic and readily computable given a set of input parameters. This means that a radio designer can choose parameters that result in the least interference with the signal of interest before committing to a specific configuration.\n",
    "\n",
    "## 2. Effects of Non-Linearity <a class=\"anchor\" id=\"non_lin\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ea202fa-d90b-46a6-806f-e140a40ab281",
   "metadata": {},
   "source": [
    "A linear amplifier will provide exactly linear gain across its entire operating range, meaning that the output of a device will be exactly proportionate to the input signal across all amplitudes. The output of a linear amplifier with a constant gain coefficient, $a$, can be expressed by\n",
    "\n",
    "$$y = ax.$$\n",
    "\n",
    "The output of a non-linear device on the other hand will not be proportionate for all input amplitudes and can be expressed by\n",
    "\n",
    "$$a_1x + a_2x^2 + a_3x^3 ... + a_nx^n.$$\n",
    "\n",
    "$a_1$, $a_2$... are not parameters a user can control, rather they are coefficients used to model the physical properties of the amplifier. The ideal scenario is that $a_2$ and above will be equal to zero. However, in practice they will have non-zero (but usually small) values.\n",
    "\n",
    "Let's look at how a non-linear response can affect the output of an amplifier.\n",
    "\n",
    "First we set the non-linear coefficients of the amplifier. In this example, the amplifier we are modeling has a 10x gain, which will be the value of $a_1$. We will only be using $a_2$ and $a_3$ as the non-linear coefficients. In practice, these would continue on indefinitely (usually with the higher order coefficient values tending towards zero). Note that the values we use for $a_2$ and $a_3$ are exaggerated for demonstration purposes, modeling a *very* non-linear amplifier, with real values likely to be much lower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f301ce8-173d-4102-9c60-0ac5a2463d41",
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = 10\n",
    "a2 = -2.2\n",
    "a3 = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd869ae4-9a3f-4eac-b172-abcf86f9ee96",
   "metadata": {},
   "source": [
    "We can now set up the inputs and outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "864e0b54-5fce-40d7-8c74-aa7dbda321cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x = np.arange(-1,1,0.001)\n",
    "y1 = a1*x\n",
    "y2 = a1*x + a2*x**2 + a3*x**3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df4b10b-7858-4fb9-9383-4a2123ad2564",
   "metadata": {},
   "source": [
    "Running the next cell will plot the input and output amplitudes for both linear and non-linear responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c94131b6-9c67-4e2a-be7b-5aa335bcabaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.plot(x, y1, label='Linear')\n",
    "plt.plot(x, y2, label='Non-Linear')\n",
    "\n",
    "plt.xlabel('Input signal amplitude')\n",
    "plt.ylabel('Output signal amplitude')\n",
    "plt.xlim(-1,1)\n",
    "plt.ylim(-15,15)\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74dd81fd-5587-4e2f-9d5d-2bddbb4755ae",
   "metadata": {},
   "source": [
    "As we can see from the plot, the output of the linear response provides exactly linear gain, whereas the non-linear response diverges as the input signal amplitude increases. This divergence will cause spurs to occur, with frequency components harmonically related to the input signal.\n",
    "\n",
    "We can test this out further by simulating a non-linear device with a sinusoid input signal.\n",
    "\n",
    "First we need to define the parameters and create the input signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ce2dcdf-8dd8-436b-9548-8829284dbaeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 10e3 # sample rate\n",
    "T = np.arange(0,1,1/fs) \n",
    "\n",
    "Nfft = 1024 # number of FFT points.\n",
    "fin = 100*(fs/Nfft) # input frequency\n",
    "\n",
    "sin_x = np.sin(2*np.pi*fin*T)\n",
    "\n",
    "print(\"fs: {:.2f} Hz\".format(fs))\n",
    "print(\"f_in: {:.2f} Hz\".format(fin))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f60513c-1ed9-4498-85b3-5a42b478a9ed",
   "metadata": {},
   "source": [
    "The input signal is a sine wave with a frequency of 976.56 Hz and a sample rate of 10 kHz. We will be using these parameters throughout the rest of this notebook.\n",
    "\n",
    "We then create two output waveforms, one with a linear response, the other with a non-linear response. Note we are using the same coefficients we defined earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dc82ff1-61bc-46dc-ae60-a758847acd3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sin_y_lin = a1*sin_x\n",
    "sin_y_nonlin = a1*sin_x + a2*sin_x**2 + a3*sin_x**3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1d62adf-4c58-42fd-83aa-8c0ba07b3959",
   "metadata": {},
   "source": [
    "We can then plot the waveforms and compare them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9fb1430-2d4a-40af-bbd0-3807e90f141a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.subplot(3,1,1)\n",
    "plt.plot(T, sin_x, label='Input Sine Wave', color='tab:blue')\n",
    "plt.xlim(0,0.01)\n",
    "plt.ylim(-15,15)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend(loc='upper right')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplot(3,1,2)\n",
    "plt.plot(T, sin_y_lin, label='Output (Linear)', color='tab:orange')\n",
    "plt.xlim(0,0.01)\n",
    "plt.ylim(-15,15)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend(loc='upper right')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.subplot(3,1,3)\n",
    "plt.plot(T, sin_y_nonlin, label='Output (Non-Linear)', color='tab:red')\n",
    "plt.xlim(0,0.01)\n",
    "plt.ylim(-15,15)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend(loc='upper right')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21f98552-a49b-40cd-a9de-1fe942492838",
   "metadata": {},
   "source": [
    "As expected, the waveform with the linear response provides linear gain to the input signal. The non-linear response, however, has changed the signal slightly. \n",
    "\n",
    "While the difference between the linear and non-linear responses may be difficult to spot in the time domain, in the frequency domain it is a different matter.\n",
    "\n",
    "In the cell below, we perform an FFT on both output signals to compare them in the frequency domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17ebcfc9-cc50-489b-8abc-5276966f7bd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_FFT = np.linspace(-fs/2, fs/2, Nfft, endpoint=False)\n",
    "\n",
    "Y_nonlin = abs(np.fft.fftshift(np.fft.fft(sin_y_nonlin, Nfft)))\n",
    "Y_nonlin_norm = Y_nonlin/Y_nonlin.max()\n",
    "Y_lin =  abs(np.fft.fftshift(np.fft.fft(sin_y_lin, Nfft)))\n",
    "Y_lin_norm = Y_lin/Y_lin.max()\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.subplot(211)\n",
    "plt.plot(F_FFT, 20*np.log10(Y_lin_norm), label='Linear', color='tab:blue')\n",
    "plt.xlim(0,fs/2)\n",
    "plt.ylabel(\"Magnitude (dBc)\")\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.grid(True)\n",
    "plt.legend(loc='upper right')\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.plot(F_FFT, 20*np.log10(Y_nonlin_norm), label='Non-Linear', color='tab:orange')\n",
    "plt.xlim(0,fs/2)\n",
    "plt.ylabel(\"Magnitude (dBc)\")\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.grid(True)\n",
    "plt.legend(loc='upper right')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "942698a3-4cb0-4889-8c3d-78fc25f068a5",
   "metadata": {},
   "source": [
    "As we can see, the simulated non-linear device has frequency components at 976.56 Hz (the fundamental), 1.95 kHz (the second harmonic), 2.93 kHz (the third harmonic), as well as a component at 0 Hz.\n",
    "\n",
    "Manufacturer data sheets typically only include the second and third harmonics as they tend to be the most significant. However, higher order harmonics will exist to some extent in all non-linear devices, albeit usually at lower amplitudes.\n",
    "\n",
    "We can easily calculate the location of any harmonic using the following equation:\n",
    "\n",
    "$$HD_n = f_{in}n,$$\n",
    "\n",
    "where $f_{in}$ is the input frequency and $n$ is the order of the harmonic (2, 3, 4...).\n",
    "\n",
    "In the cell below, we use the equation above to calculate the second and third harmonics and plot them on a frequency plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b35df81a-781f-4105-a246-e6bda7202ed3",
   "metadata": {},
   "outputs": [],
   "source": [
    "hd_ord = [2,3]\n",
    "\n",
    "hd = [fin*n for n in hd_ord]\n",
    "print(\"f_in: {:.2f} Hz\".format(fin))\n",
    "print(\"Harmonics: {} (Hz)\".format(hd))\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.vlines(fin, ymin=0, ymax=1, label='Input Signal', color='tab:blue')\n",
    "plt.vlines(hd, ymin=0, ymax=0.5, label='Harmonics', color='tab:orange')\n",
    "plt.vlines(fs/2, ymin=0, ymax=1.1, label='Nyquist Rate', color='black', linestyles='dashed')\n",
    "plt.vlines(fs, ymin=0, ymax=1.1, label='Sample Rate', color='gray', linestyles='dashed')\n",
    "\n",
    "plt.ylim(0,1.1)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Magnitude (normalised)\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c5f3cba-d94b-4ee6-bf7a-efe8b40ca4b2",
   "metadata": {},
   "source": [
    "As expected, both harmonics are at the same frequencies as the simulated device. Note that in the plot above the \"magnitude\" of the spurs is for demonstrative purposes only. The actual amplitude can only be found by direct measurement."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ceab1663-9ff2-4f03-b0e7-0339870db86f",
   "metadata": {},
   "source": [
    "## 3. ADC Interleaving Spurs <a class=\"anchor\" id=\"il_spurs\"></a>\n",
    "\n",
    "As discussed in Chapter 12, spurs can occur in interleaving ADCs due to mismatches between the sub-ADCs. The most prominent are DC offset, gain, and phase mismatches. \n",
    "\n",
    "As we did with a non-linear device before, we can simulate an interleaving ADC and see how these mismatches can introduce spurs to the output signal.\n",
    "\n",
    "First we need to create a \"continuous\" input signal the ADC can sample. We do this by using a sample rate much higher than the original signal we used earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82cb7b0e-5412-485f-b1f3-5f7be6fa6a5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs_cont = fs*10\n",
    "t_cont = np.arange(0,1,1/fs_cont)\n",
    "\n",
    "sin_y_cont = np.sin(2*np.pi*fin*t_cont)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a91835d-7a67-42ed-b07b-7f251e6b0f5e",
   "metadata": {},
   "source": [
    "Now we have a \"continuous\" signal, we can set some parameters for the interleaving ADCs.\n",
    "\n",
    "For this example, we will be using an interleaving factor, $M$, of 2. Additionally, we will set the sample rate of the sub-ADCs to $f_s/M = 10^3/2 = 5 kHz$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4b7c919-01cc-44f3-b233-f708a7b9aa09",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 2 # interleaving factor\n",
    "fs_sub_adc = fs/M # sub ADC sample rate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee2a4560-c6ab-443c-8f59-247ef2202b19",
   "metadata": {},
   "source": [
    "Since each sub-ADC needs to sample at different intervals, we need to calculate their individual sample points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ce80e9-ab7c-4ab3-b028-0e255f2bd229",
   "metadata": {},
   "outputs": [],
   "source": [
    "adc1_t_offset = round((fs_cont/fs_sub_adc)*((1-1)/M))\n",
    "adc2_t_offset = round((fs_cont/fs_sub_adc)*((2-1)/M))  \n",
    "\n",
    "t_sub_adc1 = np.arange(adc1_t_offset,len(sin_y_cont),fs_cont/fs_sub_adc, dtype=np.int16)\n",
    "t_sub_adc2 = np.arange(adc2_t_offset,len(sin_y_cont),fs_cont/fs_sub_adc, dtype=np.int16)\n",
    "\n",
    "print(\"ADC1 Sample Points: {}\".format(t_sub_adc1[0:5]))\n",
    "print(\"ADC2 Sample Points: {}\".format(t_sub_adc2[0:5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ed6e4b-e47b-4ac5-b234-9e530739eebd",
   "metadata": {},
   "source": [
    "We can see from the first 5 sample points of each ADC that they alternate.\n",
    "\n",
    "Next we define the mismatch offsets between the sub-ADCs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "908fc5fe-565a-40cc-8a7e-c86ff5ba54e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "dc_offset = 0.2\n",
    "gain_offset = 0.4\n",
    "phase_offset = 0.3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c58c003-701c-4e56-8fa7-e42d2ba7367d",
   "metadata": {},
   "source": [
    "We then sample the continuous signal at the sample points for each sub-ADC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b16c58f5-59d4-454c-bad3-0c78e1470675",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_adc1_out = sin_y_cont[t_sub_adc1]\n",
    "sub_adc2_out = sin_y_cont[t_sub_adc2 + int(np.round(phase_offset * adc2_t_offset))] * (1-gain_offset) + dc_offset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3622687-e183-4104-a61d-61502b6c8cb8",
   "metadata": {},
   "source": [
    "Finally, we interleave the outputs from each sub-ADC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a79ee696-5607-4c20-90d4-306b271ff3d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "adc_out = [None]*(len(sub_adc1_out) + len(sub_adc2_out))  \n",
    "adc_out[::M] = sub_adc1_out\n",
    "adc_out[1::M] = sub_adc2_out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d26cbb0-a444-40e4-9407-915e322773a1",
   "metadata": {},
   "source": [
    "We can now take the FFT of the interleaved output signal and view it in the frequency domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5fc820e-93e8-4b7b-8131-1bba2833c067",
   "metadata": {},
   "outputs": [],
   "source": [
    "ADC_OUT = abs(np.fft.fftshift(np.fft.fft(adc_out, Nfft)))\n",
    "ADC_OUT_NORM = ADC_OUT/ADC_OUT.max()\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.plot(F_FFT,20*np.log10(ADC_OUT_NORM))\n",
    "\n",
    "plt.xlim(-10,fs/2 + 10)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Magnitude (dBc)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "224ba294-91a0-4b9d-984f-ae09ac4076db",
   "metadata": {},
   "source": [
    "As shown in the above plot, the interleaved ADC we have simulated has frequency components at 976.56 Hz (the fundamental), 4.02 kHz (gain/time mismatch), and 0 Hz (DC offset interleaving spur). There should also be a DC offset spur at exactly $f_s/2$. However, due to the way we are plotting the frequency domain, this spur also appears at 0 Hz.\n",
    "\n",
    "In the cell below, we formalise this simulation of the interleaved ADC into a generic function, which we will use later on in this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "208f9064-9f6e-4e71-81a4-c62cabffd80c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def interleave_adc(M, fs_cont, fs_sub_adc, sig, dc_offset, gain_offset, phase_offset):\n",
    "    \n",
    "    if (len(dc_offset) and len(gain_offset) and len(phase_offset)) != M:\n",
    "        raise RuntimeError(\"Offset values must be a list with length equal to the interleaving factor.\")\n",
    "\n",
    "    adc_t_offset = round(fs_cont/(fs_sub_adc*M))    \n",
    "    sub_adc = [[]]*M \n",
    "    adc_len = 0\n",
    "    \n",
    "    for j in range(M):\n",
    "        \n",
    "        t = np.arange(j*adc_t_offset,len(sig),fs_cont/fs_sub_adc)\n",
    "\n",
    "        for n in t:\n",
    "            sub_adc[j] = sub_adc[j] + [sig[int(n) + (int(phase_offset[j] * adc_t_offset))] * (1-gain_offset[j]) + dc_offset[j]]\n",
    "        adc_len += len(sub_adc[j])\n",
    "\n",
    "    adc = [None]*adc_len\n",
    "    \n",
    "    for i in range(M):\n",
    "        adc[i::M] = sub_adc[i]\n",
    "\n",
    "    return np.array(adc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0015850-8614-43e7-9170-08820ad8fec4",
   "metadata": {},
   "source": [
    "Similar to harmonics, interleaving spurs are relatively easy to calculate. \n",
    "\n",
    "DC offset spurs (OIS) can be calculated by:\n",
    "\n",
    "$$\\frac{k}{M}f_s,$$\n",
    "\n",
    "where $f_s$ is the sample rate, $M$ is the total number of interleaving ADCs, and $k$ is an integer $0,1,2..M-1$.\n",
    "\n",
    "Gain and time interleaving spurs have a similar relationship but are also dependent on the frequency of the input signal.\n",
    "\n",
    "$$\\frac{k}{M}f_s \\pm f_{in}.$$\n",
    "\n",
    "Additionally, any harmonics present in the signal, caused by a non-linear response, will also be affected by the interleaving ADCs and introduce spurs. These can be calculated by\n",
    "\n",
    "$$\\frac{k}{M}f_s \\pm HD_{n},$$\n",
    "\n",
    "where $HD_n$ is the $n^{th}$ order harmonic.\n",
    "\n",
    "In the following three sections we use these equations to calculate the spurs for a 4-interleaved ADC.\n",
    "\n",
    "First, we need to change the interleaving factor, $M$, to 4, and update the sample rate of the sub-ADCs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "272f71da-db8e-4c81-bacd-afec150ae164",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 4 # change interleaving factor\n",
    "fs_sub_adc = fs/M # recalculate sub-ADC sample rate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "258d1769-b4c5-45d2-9425-147e39cd1f38",
   "metadata": {},
   "source": [
    "### 3.1. DC Offset Interleaving Spurs <a class=\"anchor\" id=\"dc_offset\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4b87441-17a8-481c-9aed-21e2e00ef949",
   "metadata": {},
   "outputs": [],
   "source": [
    "ois = [k/M * fs for k in range(M)]\n",
    "print(\"fin: {:.2f} Hz\".format(fin))\n",
    "print(\"OIS: {} (Hz)\".format(ois))\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.vlines(fs/2, ymin=0, ymax=1.1, label='Nyquist Rate', color='black', linestyles='dashed')\n",
    "plt.vlines(fs, ymin=0, ymax=1.1, label='Sample Rate', color='black', linestyles='dashdot')\n",
    "plt.vlines(fin, ymin=0, ymax=1, label='Input Signal', color='tab:blue')\n",
    "plt.vlines(ois, ymin=0, ymax=0.4, label='OIS', color='tab:orange')\n",
    "\n",
    "plt.xlim(-100,fs+100)\n",
    "plt.ylim(0,1.1)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b60fa7f6-8726-4d6b-9df9-29dc70bc12db",
   "metadata": {},
   "source": [
    "### 3.2. Gain/Time Interleaving Spurs <a class=\"anchor\" id=\"gain_time\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e872ea6-f4a4-4b93-b08b-fa0747871243",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtis_p = [k*fs/M + fin for k in range(1,M)]\n",
    "gtis_n = [k*fs/M - fin for k in range(1,M)]\n",
    "\n",
    "gtis = gtis_p + gtis_n\n",
    "print(\"fin: {:.2f} Hz\".format(fin))\n",
    "print(\"GTIS: {} (Hz)\".format(gtis))\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.vlines(fs/2, ymin=0, ymax=1.1, label='Nyquist Rate', color='black', linestyles='dashed')\n",
    "plt.vlines(fs, ymin=0, ymax=1.1, label='Sample Rate', color='black', linestyles='dashdot')\n",
    "plt.vlines(fin, ymin=0, ymax=1, label='Input Signal', color='tab:blue')\n",
    "plt.vlines(gtis, ymin=0, ymax=0.3, label='GTIS', color='tab:orange')\n",
    "\n",
    "plt.xlim(-100,fs+100)\n",
    "plt.ylim(0,1.1)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.legend(loc='upper right')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b89f8034-0677-4d95-90da-9b9803953c9b",
   "metadata": {},
   "source": [
    "### 3.3. Harmonic Interleaving Spurs <a class=\"anchor\" id=\"harmonic\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27c45589-78d0-4727-925d-ca3077d7c20c",
   "metadata": {},
   "outputs": [],
   "source": [
    "hd2_il_p = [spur + hd[0] for spur in ois]\n",
    "hd2_il_n = [spur - hd[0] for spur in ois]\n",
    "\n",
    "hd3_il_p = [spur + hd[1] for spur in ois]\n",
    "hd3_il_n = [spur - hd[1] for spur in ois]\n",
    "\n",
    "hd_il = hd2_il_p + hd2_il_n + hd3_il_p + hd3_il_n\n",
    "hd_il = [abs(spur) for spur in hd_il]\n",
    "\n",
    "print(\"fin: {:.2f} Hz\".format(fin))\n",
    "print(\"HD+IL: {} (Hz)\".format(hd_il))\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.vlines(fs/2, ymin=0, ymax=1.1, label='Nyquist Rate', color='black', linestyles='dashed')\n",
    "plt.vlines(fs, ymin=0, ymax=1.1, label='Sample Rate', color='black', linestyles='dashdot')\n",
    "plt.vlines(fin, ymin=0, ymax=1, label='Input Signal', color='tab:blue')\n",
    "plt.vlines(hd_il, ymin=0, ymax=0.2, label='Harmonic + IL Spurs', color='tab:orange')\n",
    "\n",
    "plt.xlim(-100,fs+100)\n",
    "plt.ylim(0,1.1)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.legend(loc='upper right')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f795e3d-be24-420d-bbd5-b160300f5089",
   "metadata": {},
   "source": [
    "As previously noted, the amplitudes of the spurs in these plots are for demonstrative purposes only in order to make it easier to differentiate between them. The actual amplitude of the spurs can only be determined through direct measurement.\n",
    "\n",
    "## 4. Aliasing <a class=\"anchor\" id=\"alias\"></a>\n",
    "\n",
    "ADC-related spurs with frequency components higher than the Nyquist rate ($f_s/2$) will fold back into the first Nyquist zone, increasing the likelihood that a spur will interfere with the signal of interest. A spur located in the first Nyquist zone will not be aliased and so its frequency content will not be changed, while a spur located in the second Nyquist zone will alias at a frequency of $f_s - f_{in}$. The spectrum between $0$ to $f_s$ is identical to that between $f_s$ to $2f_s$ and $2_fs$ to $3f_s$ and so on. We can use this information to easily calculate the location of any spur after aliasing.\n",
    "\n",
    "Below, we define a function to determine the frequency of a spur after aliasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e06351d-438b-43f9-9e15-d8ce41eba5e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def alias(f, fs):\n",
    "    f = abs(f) # any negative frequencies will also appear as positive frequencies\n",
    "    f_alias = f % fs\n",
    "    if f_alias > fs/2:\n",
    "        f_alias = fs - f_alias\n",
    "    return f_alias"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8ae9f35-bf91-4491-b15c-d7b5813f6669",
   "metadata": {},
   "source": [
    "We can then use this function to calculate the location of all the spurs we have calculated so far."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e26970b4-f5ac-48b3-bc44-c6897bb133ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "hd_alias = [alias(f, fs) for f in hd]\n",
    "ois_alias = [alias(f, fs) for f in ois]\n",
    "gtis_alias = [alias(f, fs) for f in gtis]\n",
    "hd_il_alias = [alias(f, fs) for f in hd_il]\n",
    "\n",
    "hd_alias = list(set(hd_alias))\n",
    "ois_alias = list(set(ois_alias))\n",
    "gtis_alias = list(set(gtis_alias))\n",
    "hd_il_alias = list(set(hd_il_alias))\n",
    "\n",
    "print(\"fin: {:.2f} Hz\".format(fin))\n",
    "print(\"HD Alias: {} (Hz)\".format(hd_alias))\n",
    "print(\"OIS Alias: {} (Hz)\".format(ois_alias))\n",
    "print(\"GTIS Alias: {} (Hz)\".format(gtis_alias))\n",
    "print(\"HD+IL Alias: {} (Hz)\".format(hd_il_alias))\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.vlines(fs/2, ymin=0, ymax=1.1, label='Nyquist Rate', color='black', linestyles='dashed')\n",
    "plt.vlines(fs, ymin=0, ymax=1.1, label='Sample Rate', color='black', linestyles='dashdot')\n",
    "\n",
    "plt.vlines(fin, ymin=0, ymax=1, label='Input Signal', color='tab:blue')\n",
    "plt.vlines(hd_alias, ymin=0, ymax=0.5, label='Harmonics', color='tab:orange')\n",
    "plt.vlines(ois_alias, ymin=0, ymax=0.4, label='OIS', color='tab:green')\n",
    "plt.vlines(gtis_alias, ymin=0, ymax=0.3, label='GTIS', color='tab:cyan')\n",
    "plt.vlines(hd_il_alias, ymin=0, ymax=0.2, label='Harmonic + IL Spurs', color='tab:pink')\n",
    "\n",
    "plt.xlim(-100,fs+100)\n",
    "plt.ylim(0,1.1)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.legend(loc='upper right')\n",
    "plt.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89194982-97b5-44d1-808a-ff29d2d820ea",
   "metadata": {},
   "source": [
    "As we can see from this plot, the spectrum is now very busy, but there is ample space between the spurs closest to the input signal to filter them out."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdd57f76-46cf-4d07-8f9a-1b44dfd72df1",
   "metadata": {},
   "source": [
    "## 5. Simulating a Non-Linear Interleaved ADC <a class=\"anchor\" id=\"sim\"></a>\n",
    "\n",
    "Now that we understand how and why spurs occur, we can simulate an entire non-linear interleaved ADC and check the results against the spurs we calculated in the previous sections.\n",
    "\n",
    "First we need to define the sub-ADC offsets. Remember, the values we are defining here are exaggerated for demonstrative purposes. In practice these values will typically be much lower, and the amplitude of spurs can only be accurately determined by direct measurement.\n",
    "\n",
    "Since we have 4 sub-ADCs now we need to define a mismatch for each one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d920bed-4304-463f-a5bb-955b3c86c5e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "dc_offset = [0, 0.01, 0.05, 0.02]\n",
    "gain_offset = [0, 0.01, 0.02, 0.03]\n",
    "phase_offset = [0, 0.1, 0.2, 0.01]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95c905db-04fc-4ace-aceb-0764368b1dc9",
   "metadata": {},
   "source": [
    "Then we can calculate the output of the simulated interleaved ADC using the `interleave_adc()` function we defined earlier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a4654ef-c10d-47a9-b841-2ed2748a95f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "adc_out = interleave_adc(M, fs_cont, fs_sub_adc, sin_y_cont, dc_offset, gain_offset, phase_offset)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5af26e8d-c774-49e6-8fec-e20f759e5c1d",
   "metadata": {},
   "source": [
    "We can then simulate the response of the non-linear ADC.\n",
    "\n",
    "First we need to redefine the coefficients. Note that we are using a value of 1 for `a1` as the ADC does not amplify the signal and we should expect an equivalent amplitude between input and output for a linear device."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a10e4fb8-98c3-4f94-99f6-6193bffff67b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = 1\n",
    "a2 = 0.05\n",
    "a3 = 0.02"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f785a6f-7a76-4c35-9511-cf7480b1914f",
   "metadata": {},
   "source": [
    "We can then use the equation we defined earlier to simulate a non-linear response from the ADC. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ea2847e-c344-4748-b768-56d6837b4d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "adc_out_nonlin = a1*adc_out + a2*adc_out**2 + a3*adc_out**3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a735a42-7034-4423-b3cc-70035185fdd6",
   "metadata": {},
   "source": [
    "> It should be noted that this is not a true simulation of a non-linear ADC. The gain differentials between the input and output of an ADC occur from slightly different, but related, phenomena than in an amplifier. However, the end result is the same and for the purposes of this notebook we can calculate the non-linear response using the equation above. \n",
    "\n",
    "We can also add a little noise to the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9d33327-250c-43f0-8f77-7c9ce205dfff",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise = np.random.normal(0,1,len(adc_out_nonlin))\n",
    "adc_out_nonlin += noise*0.0001"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1424bc37-56d7-4a07-ab16-7c48b6f542cd",
   "metadata": {},
   "source": [
    "Finally, we can take the FFT of the output signal and overlay the location of the aliased spurs we calculated a few cells above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea65b83f-3d73-4229-8248-21f9e6d2b98e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ADC_OUT = abs(np.fft.fftshift(np.fft.fft(adc_out_nonlin, Nfft)))\n",
    "ADC_OUT_NORM = ADC_OUT/ADC_OUT.max()\n",
    "\n",
    "alpha = 0.8\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.vlines(ois_alias, ymin=-300, ymax=0, label='OIS', color='tab:green', alpha=alpha, linestyles='dashed')\n",
    "plt.vlines(gtis_alias, ymin=-300, ymax=0, label='GTIS', color='tab:cyan', alpha=alpha, linestyles='dashed')\n",
    "plt.vlines(hd_il_alias, ymin=-300, ymax=0, label='Harmonic + IL Spurs', color='tab:pink', alpha=alpha, linestyles='dashed')\n",
    "plt.vlines(hd_alias, ymin=-300, ymax=0, label='Harmonics', color='tab:orange', alpha=alpha, linestyles='dashed')\n",
    "plt.plot(F_FFT,20*np.log10(ADC_OUT_NORM), color='tab:blue', label='ADC Output')\n",
    "\n",
    "plt.xlim(-10,fs/2 + 10)\n",
    "plt.ylim(-320,10)\n",
    "plt.ylim(-115,10)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Magnitude (dBc)\")\n",
    "plt.legend(loc=\"lower right\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5892c020-f6ad-48f8-8b0d-6d6e853b703e",
   "metadata": {},
   "source": [
    "As we can see from the plot, the location of the spurs are exactly where we had calculated them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94b24a84-7a4f-4a78-a0af-0ce2df4b40ef",
   "metadata": {},
   "source": [
    "## 6. Filter Design <a class=\"anchor\" id=\"filt_design\"></a>\n",
    "\n",
    "Now that we have our simulated output we can design a filter to suppress the surrounding spurs. First we need to determine the transition band for the filter. \n",
    "\n",
    "We can see from the above plot that the closest spurs to the input signal are the GTIS and harmonic interleaving spurs. We can print out these values to determine exactly how close they are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "335feb71-b5b9-445f-9493-75abcd2f2c4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtis_alias.sort()\n",
    "hd_il_alias.sort()\n",
    "\n",
    "print(\"f_in: {} (Hz)\".format(fin))\n",
    "print(\"HD_IL: {} (Hz)\".format(hd_il_alias))\n",
    "print(\"GTIS: {} (Hz)\".format(gtis_alias))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb4b1265-2638-4b24-882d-3b225b8206ff",
   "metadata": {},
   "source": [
    "The closest harmonic interleaving spurs is at around 546 Hz, while the closest GTIS spur is around 1.5 kHz. A transition band of around 400 Hz should be more than sufficient.\n",
    "\n",
    "Since there are spurs on either side of the input signal we will require a bandpass filter. \n",
    "\n",
    "Below we design a filter using the `firwin2` function from the SciPy library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7bcece5-b3a9-4ed9-a5ae-cfaa9c641a49",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import firwin2\n",
    "\n",
    "num_taps = 51\n",
    "bpf = firwin2(num_taps, [0,fin-400,fin,fin+400,fs/2], [0,0,1,0,0], fs=fs, window='blackman')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "444446b8-8862-4df0-b590-dea34da8ab5c",
   "metadata": {},
   "source": [
    "We can then view the response of filter in the frequency domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19cdf6a3-8398-4c12-81f9-8818c653100a",
   "metadata": {},
   "outputs": [],
   "source": [
    "BPF = np.fft.fftshift(np.fft.fft(bpf, Nfft))\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.plot(F_FFT,20*np.log10(abs(BPF)))\n",
    "\n",
    "plt.xlim(-10,fs/2 + 10)\n",
    "plt.ylim(-150,10)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Magnitude (dB)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e965b85b-9396-4ecf-963f-ca76019fe40c",
   "metadata": {},
   "source": [
    "Finally we can apply the filter to the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40f6ad01-2ac7-43bf-832e-d2b5c878d6b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "adc_out_filt = np.convolve(adc_out_nonlin, bpf, 'same')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d22d734-b066-42e9-9992-56cb46119f5e",
   "metadata": {},
   "source": [
    "We then take the FFT and plot the results in the frequency domain. We will overlay the original signal to see how well the filter has worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43566ec3-2eca-436c-9fc5-d72b3e72f57e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ADC_OUT_FILT = abs(np.fft.fftshift(np.fft.fft(adc_out_filt, Nfft)))\n",
    "ADC_OUT_FILT_NORM = ADC_OUT_FILT/ADC_OUT_FILT.max()\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.plot(F_FFT,20*np.log10(ADC_OUT_NORM), color='tab:orange', label='Unfiltered ADC output')\n",
    "plt.plot(F_FFT,20*np.log10(ADC_OUT_FILT_NORM), color='tab:blue', label='Filtered ADC output')\n",
    "\n",
    "plt.xlim(-10,fs/2 + 10)\n",
    "plt.ylim(-115,5)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Magnitude (dBc)\")\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b7c27dd-94da-41aa-b967-f8dc9ad33cc6",
   "metadata": {},
   "source": [
    "We can see from the plot that the filter has sufficiently removed the spurs from the signal, leaving around -50 dBc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "655f9a59-9a07-4600-bc73-3e9140bed48a",
   "metadata": {},
   "source": [
    "## 7. Conclusion <a class=\"anchor\" id=\"conclusion\"></a>\n",
    "\n",
    "In this notebook we dove a little deeper into how and why spurs occur within data converters and how to calculate where they will appear on the spectrum. We showed how to simulate the effects of non-linearity and interleaving, and compared the simulated results to the calculated ones.\n",
    "\n",
    "In the next notebook, we present a frequency planning tool that uses the concepts introduced in this notebook to allow the user to quickly and efficiently calculate spurs given arbitrary input parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc1a4073-2775-47e4-8e51-a3377b51abb6",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[⬅️ Previous Notebook](../notebook_E/03_complex_qam.ipynb) || [Next Notebook 🚀](02_rfsoc_frequency_planner.ipynb)\n",
    "\n",
    "Copyright © 2023 Strathclyde Academic Media\n",
    "\n",
    "---\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
